\defmodule{MultinomialDist}

Implements the abstract class \class{DiscreteDistributionIntMulti} for the
{\em multinomial} distribution with parameters $n$ and
($p_1$, \ldots,$p_d$).
The probability mass function is \cite{tJOH69a}
\begin{htmlonly}
\eq
   P[X = (x_1,...,x_d)] = {n!} \prod_{i=1}^{d}p_i^{x_i}/x_i!,
\endeq
\end{htmlonly}
\begin{latexonly}
\eq
   P[X = (x_1,\ldots,x_d)] = {n!} \prod_{i=1}^{d}\frac{p_i^{x_i}}{x_i!},
\eqlabel{eq:fMultinomial}
\endeq
\end{latexonly}
where $\sum_{i=1}^{d} x_i = n$ and $\sum_{i=1}^{d} p_i = 1$.

\bigskip\hrule

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{code}
\begin{hide}
/*
 * Class:        MultinomialDist
 * Description:  multinomial distribution
 * Environment:  Java
 * Software:     SSJ
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author
 * @since

 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.

 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.

 * A copy of the GNU General Public License is available at
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */
\end{hide}
package umontreal.iro.lecuyer.probdistmulti;
\begin{hide}
import umontreal.iro.lecuyer.util.Num;
\end{hide}

public class MultinomialDist extends DiscreteDistributionIntMulti \begin{hide} {
   protected int n;
   protected double p[];

\end{hide}
\end{code}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection* {Constructors}

\begin{code}

   public MultinomialDist (int n, double p[]) \begin{hide} {
      setParams (n, p);
   }\end{hide}
\end{code}
\begin{tabb}
   Creates a \texttt{MultinomialDist} object with parameters $n$ and
   ($p_1$,\ldots,$p_d$) such that $\sum_{i=1}^{d} p_i = 1$. We have
   $p_i = $ \texttt{p[i-1]}.
\end{tabb}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection* {Methods}

\begin{code}\begin{hide}
   public double prob (int x[]) {
      return prob_ (n, p, x);
   }

   public double cdf (int x[]) {
      return cdf_ (n, p, x);
   }

   public double[] getMean() {
      return getMean_ (n, p);
   }

   public double[][] getCovariance() {
      return getCovariance_ (n, p);
   }

   public double[][] getCorrelation () {
      return getCorrelation_ (n, p);
   }

   private static void verifParam(int n, double p[]) {
      if (n <= 0)
         throw new IllegalArgumentException ("n <= 0");

      double sumPi = 0.0;
      for (int i = 0; i < p.length; i++) {
         if ((p[i] < 0) || (p[i] > 1))
            throw new IllegalArgumentException("p is not a probability vector");
         sumPi += p[i];
      }

      if (sumPi != 1.0)
         throw new IllegalArgumentException ("p is not a probability vector");
   }

   private static double prob_ (int n, double p[], int x[]) {
      if (x.length != p.length)
         throw new IllegalArgumentException ("x and p must have the same dimension");

      double sumXFact = 0.0;
      int sumX = 0;
      double sumPX = 0.0;

      for (int i = 0; i < p.length; i++) {
         sumX += x[i];
         sumXFact += Num.lnFactorial (x[i]);
         sumPX += (x[i] * Math.log (p[i]));
      }

      if (sumX != n)
         return 0.0;
      else {
         return Math.exp (Num.lnFactorial (n) - sumXFact + sumPX);
      }
   }
\end{hide}

   public static double prob (int n, double p[], int x[])\begin{hide} {
      verifParam (n, p);
      return prob_ (n, p, x);
   }\end{hide}
\end{code}
\begin{tabb}
   Computes the probability mass function (\ref{eq:fMultinomial})
   of the multinomial distribution with parameters $n$ and
   ($p_1$,\ldots,$p_d$) evaluated at $x$.
\end{tabb}
\begin{code}\begin{hide}

   private static double cdf_ (int n, double p[], int x[]) {
      boolean end = false;
      double sum = 0.0;
      int j;

      if (x.length != p.length)
         throw new IllegalArgumentException ("x and p must have the same dimension");

      int is[] = new int[x.length];
      for (int i = 0; i < is.length; i++)
         is[i] = 0;

      sum = 0.0;
      while (! end) {
         sum += prob (n, p, is);
         is[0]++;

         if (is[0] > x[0]) {
            is[0] = 0;
            j = 1;
            while (j < x.length && is[j] == x[j])
               is[j++] = 0;

            if (j == x.length)
               end = true;
            else
               is[j]++;
         }
      }

      return sum;
   }\end{hide}

   public static double cdf (int n, double p[], int x[])\begin{hide} {
      verifParam (n, p);
      return cdf_ (n, p, x);
   }\end{hide}
\end{code}
\begin{tabb}
   Computes the function $F$ of the multinomial distribution with
   parameters $n$ and ($p_1$,\ldots,$p_d$) evaluated at $x$.
\end{tabb}
\begin{code}\begin{hide}

   private static double[] getMean_ (int n, double[] p) {
      double mean[] = new double[p.length];

      for (int i = 0; i < p.length; i++)
         mean[i] = n * p[i];

      return mean;
   }\end{hide}

   public static double[] getMean (int n, double[] p)\begin{hide} {
      verifParam (n, p);

      return getMean_ (n, p);
   }\end{hide}
\end{code}
\begin{tabb}
   Computes the mean $E[X_i] = np_i$ of the multinomial distribution
   with parameters $n$ and ($p_1$,\ldots,$p_d$).
\end{tabb}
\begin{code}\begin{hide}

   private static double[][] getCovariance_ (int n, double[] p) {
      double cov[][] = new double[p.length][p.length];

      for (int i = 0; i < p.length; i++) {
         for (int j = 0; j < p.length; j++)
            cov[i][j] = -n * p[i] * p[j];

         cov[i][i] = n * p[i] * (1.0 - p[i]);
      }
      return cov;
   }\end{hide}

   public static double[][] getCovariance (int n, double[] p)\begin{hide} {
      verifParam (n, p);
      return getCovariance_ (n, p);
   }\end{hide}
\end{code}
\begin{tabb}
   Computes the covariance matrix of the multinomial distribution
   with parameters $n$ and ($p_1$,\ldots,$p_d$).
\end{tabb}
\begin{code}\begin{hide}

   private static double[][] getCorrelation_ (int n, double[] p) {
      double corr[][] = new double[p.length][p.length];

      for (int i = 0; i < p.length; i++) {
         for (int j = 0; j < p.length; j++)
            corr[i][j] = -Math.sqrt(p[i] * p[j] / ((1.0 - p[i]) * (1.0 - p[j])));
         corr[i][i] = 1.0;
      }
      return corr;
   }\end{hide}

   public static double[][] getCorrelation (int n, double[] p)\begin{hide} {
      verifParam (n, p);
      return getCorrelation_ (n, p);
   }\end{hide}
\end{code}
\begin{tabb}
   Computes the correlation matrix of the multinomial distribution
   with parameters $n$ and ($p_1$,\ldots,$p_d$).
\end{tabb}
\begin{code}

   public static double[] getMLE (int x[][], int m, int d, int n)\begin{hide} {
      double parameters[] = new double[d];
      double xBar[] = new double[d];
      double N = 0.0;

      if (m <= 0)
         throw new IllegalArgumentException ("m <= 0");
      if (d <= 0)
         throw new IllegalArgumentException ("d <= 0");

      for (int i = 0; i < d; i++)
         xBar[i] = 0;

      for (int v = 0; v < m; v++)
         for (int c = 0; c < d; c++)
            xBar[c] += x[v][c];

      for (int i = 0; i < d; i++)
      {
         xBar[i] = xBar[i] / (double) n;
         N += xBar[i];
      }
      if (N != (double) n)
         throw new IllegalArgumentException("n is not correct");

      for (int i = 0; i < d; i++)
         parameters[i] = xBar[i] / (double) n;

      return parameters;
   }\end{hide}
\end{code}
\begin{tabb}
   Estimates and returns the parameters [$\hat{p_i}$,\ldots,$\hat{p_d}$] of the
   multinomial distribution using the maximum likelihood method. It uses the $m$
   observations of $d$ components in table $x[i][j]$, $i = 0, 1, \ldots, m-1$
   and $j = 0, 1, \ldots, d-1$.
   \begin{detailed}%
   The equations of the maximum likelihood are defined as %in \cite{xxxxx}:
   \begin{eqnarray*}
      \hat p_i = \frac{\bar{X_i}}{N}.
   \end{eqnarray*}
   \end{detailed}%
\end{tabb}
\begin{htmlonly}
   \param{x}{the list of observations used to evaluate parameters}
   \param{m}{the number of observations used to evaluate parameters}
   \param{d}{the dimension of each observation}
   \param{n}{the number of independant trials for each series}
   \return{returns the parameters [$\hat{p_i}$,\ldots,$\hat{p_d}$]}
\end{htmlonly}
\begin{code}

   public int getN()\begin{hide} {
      return n;
   }\end{hide}
\end{code}
\begin{tabb}
   Returns the parameter $n$ of this object.
\end{tabb}
\begin{code}

   public double[] getP()\begin{hide} {
      return p;
   }\end{hide}
\end{code}
\begin{tabb}
   Returns the parameters ($p_1$,\ldots,$p_d$) of this object.
\end{tabb}
\begin{code}

   public void setParams (int n, double p[])\begin{hide} {
      double sumP = 0.0;

      if (n <= 0)
         throw new IllegalArgumentException ("n <= 0");
      if (p.length < 2)
         throw new IllegalArgumentException ("p.length < 2");

      this.n = n;
      this.dimension = p.length;
      this.p = new double[dimension];
      for (int i = 0; i < dimension; i++)
      {
         if ((p[i] < 0) || (p[i] > 1))
            throw new IllegalArgumentException("p is not a probability vector");

         this.p[i] = p[i];
         sumP += p[i];
      }

      if (sumP != 1.0)
         throw new IllegalArgumentException ("p is not a probability vector");
   }\end{hide}
\end{code}
\begin{tabb}
   Sets the parameters $n$ and ($p_1$,\ldots,$p_d$) of this object.
\end{tabb}
\begin{code}\begin{hide}
}\end{hide}
\end{code}
